{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Team Surface Velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Members**: Grace Barcheck, Canyon Breyer, Rodrigo Gómez-Fell, Trevor Hillebrand, Ben Hills, Lynn Kaluzienski, Joseph Martin, David Polashenski"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Science Advisor**: Daniel Shapero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Special Thanks**: Ben Smith, David Shean "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Motivation\n",
    "\n",
    "**Speaker: Canyon Breyer**\n",
    "\n",
    "Previous work by Marsh and Rack (2012), and Lee and others (2012),  have demonstrated the value of using satellite altimetry as a method of calculating ice surface velocity utilizing the Geoscience Laser Altimeter System (GLAS) on board ICESat. This altimetry method has several advantages over more traditional techniques due to high pointing accuracy for geo-location and an ability to measure velocity in regions that lack visible surface features (Marsh and Rack, 2012). The method also has the added benefit of dealing with tidal fluctuations without the need for a tidal correction model. The motivation for this project was to expand the methodology outlined in Marsh and Rack (2012) to the ICE-Sat2 dataset. The smaller footprint of the ICE-Sat2 mission will likely improve the overall accuracy of velocity measurements and the nature of its precise repeat passes would provide an avenue for studying temporal variations of glacier velocities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Project Objective:\n",
    "\n",
    "**Speaker: Rodrigo Gómez-Fell**\n",
    "\n",
    "Extract surface ice velocity on polar regions from ICESat-2 along track measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Goals:\n",
    "\n",
    "- Compare the capabilities of ICESat-2 to extract surface ice velocity from ice shelves and ice streams\n",
    "- Compare ICESat GLAS methodology (along track) to ICESat-2 across track\n",
    "- Use crossovers for calculating velocities and determine how the measurements compare with simple along track and across track.\n",
    "    -Does this resolve different directions of ice flow?\n",
    "- Can a surface velocity product be extracted from ATL06, or is ATL03 the more suitable product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Study Area:\n",
    "\n",
    "When looking for a study region to test our ICESat-2 velocity derivation method, we prioritized regions that **1)** included both grounded and floating ice and **2)** had a good alignment between satellite track position and overall flow direction. We found Foundation Ice Stream, a large ice stream draining into the Filchner-Ronne Ice Shelf, to meet both of these criteria. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![FIS](FISimage.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data Selection\n",
    "\n",
    "We used the ICESat-2 Land Ice Height ATL06 product and then used the MEaSUREs Antarctic Velocity Map V2 (Rignot, 2017) for validation of our derived velocities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Method\n",
    "\n",
    "**Speaker: Ben Hills**\n",
    "\n",
    "Following Marsh and Rack (2012) we used the slope of elevation for analysis, this helped amplify differences in the ice profile between repeat measurements and also removed the influence of tidal effects. This is portrayed in the figure below.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Beardmore](beardmore.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Marsh_method](marsh2.PNG)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fig.2: From Marsh and Rack 2012. Schematic of the method used to reduce the effect of oblique surface features and ice flow which is non-parallel to ICESat tracks. Black lines indicate satellite tracks, grey ticks indicate the orientation of surface features, and  ⍺ is the feature-track angle. Bottom right profile illustrates that after adjustment there is no relative feature displacement due to cross-track separation, therefore all displacement is due to ice movement in the track direction.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Our Methods:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Cross-correlation background**\n",
    "\n",
    "**Speaker: Grace Barcheck**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Test.scipy.signal.correlate on some ATL06 data from Foundation Ice Stream (FIS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy, sys, os, pyproj, glob, re, h5py\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.signal import correlate\n",
    "from astropy.time import Time\n",
    "\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Test scipy.signal.correlate\n",
    "\n",
    "Generate test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "564af9fd03b84dc99258e25122a0ae50",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'index')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dx = 0.1\n",
    "x = np.arange(0,10,dx)\n",
    "y = np.zeros(np.shape(x))\n",
    "ix0 = 30\n",
    "ix1 = 30 + 15\n",
    "y[ix0:ix1] = 1\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y,'k')\n",
    "axs[1].set_xlabel('index')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we generate a signal to correlate the test data with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4ffaf704b8ff4fbeace23b104981c487",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'black = original, blue = shifted')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imposed_offset = int(14/dx) # 14 meters, in units of samples\n",
    "\n",
    "x_noise = np.arange(0,50,dx) # make the vector we're comparing with much longer\n",
    "y_noise = np.zeros(np.shape(x_noise))\n",
    "y_noise[ix0 + imposed_offset : ix1 + imposed_offset] = 1\n",
    "\n",
    "# uncomment the line below to add noise\n",
    "# y_noise = y_noise * np.random.random(np.shape(y_noise))\n",
    "\n",
    "fig,axs = plt.subplots(1,2)\n",
    "\n",
    "axs[0].plot(x,y,'k')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x)), y, 'k')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "axs[0].plot(x_noise,y_noise, 'b')\n",
    "axs[0].set_xlabel('distance (m)')\n",
    "axs[0].set_ylabel('value')\n",
    "axs[1].plot(np.arange(len(x_noise)), y_noise,'b')\n",
    "axs[1].set_xlabel('index')\n",
    "\n",
    "fig.suptitle('black = original, blue = shifted')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Try scipy.signal.correlate:\n",
    "\n",
    "mode ='full' returns the entire cross correlation; could be 'valid' to return only non- zero-padded part\n",
    "\n",
    "method = direct (not fft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "corr = correlate(y_noise,y, mode = 'full', method = 'direct') \n",
    "norm_val = np.sqrt(np.sum(y_noise**2)*np.sum(y**2))\n",
    "corr = corr / norm_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the dimensions of corr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "corr:  (599,)\n",
      "x:  (100,)\n",
      "x:  (500,)\n"
     ]
    }
   ],
   "source": [
    "print('corr: ', np.shape(corr))\n",
    "print('x: ', np.shape(x))\n",
    "print('x: ', np.shape(x_noise))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Look at the correlation visualized in the plots below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "84f269157c4b4f169bc305a7e75ee8d8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Shift  140  samples, or  14.0  m to line up signals')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# lagvec = np.arange(0,len(x_noise) - len(x) + 1)\n",
    "lagvec = np.arange( -(len(x) - 1), len(x_noise), 1)\n",
    "shift_vec = lagvec * dx\n",
    "\n",
    "ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "best_lag = lagvec[ix_peak]\n",
    "best_shift = shift_vec[ix_peak]\n",
    "\n",
    "fig,axs = plt.subplots(3,1)\n",
    "\n",
    "axs[0].plot(lagvec,corr)\n",
    "axs[0].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[0].set_xlabel('lag (samples)')\n",
    "axs[0].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[1].plot(shift_vec,corr)\n",
    "axs[1].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "axs[1].set_xlabel('shift (m)')\n",
    "axs[1].set_ylabel('correlation coefficient')\n",
    "\n",
    "axs[2].plot(x + best_shift, y,'k')\n",
    "axs[2].plot(x_noise, y_noise, 'b--')\n",
    "axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "\n",
    "fig.suptitle(' '.join(['Shift ', str(best_lag), ' samples, or ', str(best_shift), ' m to line up signals']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A little Background on cross-correlation..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Correlation](Corr_Coeff.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Applying our method to ATL06 data\n",
    "\n",
    "**Speaker: Ben Hills**\n",
    "\n",
    "Load repeat data:\n",
    "Import readers, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! cd ..; [ -d pointCollection ] || git clone https://www.github.com/smithB/pointCollection.git\n",
    "# sys.path.append(os.path.join(os.getcwd(), '..'))\n",
    "\n",
    "#!python3 -m pip install --user git+https://github.com/tsutterley/pointCollection.git@pip\n",
    "import pointCollection as pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "moa_datapath = '/srv/tutorial-data/land_ice_applications/'\n",
    "datapath = '/home/jovyan/shared/surface_velocity/FIS_ATL06/'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Geographic setting: Foundation Ice Stream**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/jovyan/.local/lib/python3.7/site-packages/pointCollection/__init__.py\n"
     ]
    }
   ],
   "source": [
    "print(pc.__file__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-86 -81 -81 -86 -86]\n",
      "[-55 -55 -65 -65 -55]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "29ab2392457a488ca89b1bfef8fa5b7f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-887950., -356950.,  183825.,  561825.]), 'origin': 'lower'}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Mosaic of Antarctica for Foundation Ice Stream')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spatial_extent = np.array([-65, -86, -55, -81])\n",
    "lat=spatial_extent[[1, 3, 3, 1, 1]]\n",
    "lon=spatial_extent[[2, 2, 0, 0, 2]]\n",
    "print(lat)\n",
    "print(lon)\n",
    "# project the coordinates to Antarctic polar stereographic\n",
    "xy=np.array(pyproj.Proj(3031)(lon, lat))\n",
    "# get the bounds of the projected coordinates \n",
    "XR=[np.nanmin(xy[0,:]), np.nanmax(xy[0,:])]\n",
    "YR=[np.nanmin(xy[1,:]), np.nanmax(xy[1,:])]\n",
    "MOA=pc.grid.data().from_geotif(os.path.join(moa_datapath, 'MOA','moa_2009_1km.tif'), bounds=[XR, YR])\n",
    "\n",
    "# show the mosaic:\n",
    "plt.figure()\n",
    "MOA.show(cmap='gray', clim=[14000, 17000])\n",
    "plt.plot(xy[0,:], xy[1,:])\n",
    "plt.title('Mosaic of Antarctica for Foundation Ice Stream')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Load the repeat track data\n",
    "\n",
    "ATL06 reader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def atl06_to_dict(filename, beam, field_dict=None, index=None, epsg=None):\n",
    "    \"\"\"\n",
    "        Read selected datasets from an ATL06 file\n",
    "\n",
    "        Input arguments:\n",
    "            filename: ATl06 file to read\n",
    "            beam: a string specifying which beam is to be read (ex: gt1l, gt1r, gt2l, etc)\n",
    "            field_dict: A dictinary describing the fields to be read\n",
    "                    keys give the group names to be read, \n",
    "                    entries are lists of datasets within the groups\n",
    "            index: which entries in each field to read\n",
    "            epsg: an EPSG code specifying a projection (see www.epsg.org).  Good choices are:\n",
    "                for Greenland, 3413 (polar stereographic projection, with Greenland along the Y axis)\n",
    "                for Antarctica, 3031 (polar stereographic projection, centered on the Pouth Pole)\n",
    "        Output argument:\n",
    "            D6: dictionary containing ATL06 data.  Each dataset in \n",
    "                dataset_dict has its own entry in D6.  Each dataset \n",
    "                in D6 contains a numpy array containing the \n",
    "                data\n",
    "    \"\"\"\n",
    "    if field_dict is None:\n",
    "        field_dict={None:['latitude','longitude','h_li', 'atl06_quality_summary'],\\\n",
    "                    'ground_track':['x_atc','y_atc'],\\\n",
    "                    'fit_statistics':['dh_fit_dx', 'dh_fit_dy']}\n",
    "    D={}\n",
    "    # below: file_re = regular expression, it will pull apart the regular expression to get the information from the filename\n",
    "    file_re=re.compile('ATL06_(?P<date>\\d+)_(?P<rgt>\\d\\d\\d\\d)(?P<cycle>\\d\\d)(?P<region>\\d\\d)_(?P<release>\\d\\d\\d)_(?P<version>\\d\\d).h5')\n",
    "    with h5py.File(filename,'r') as h5f:\n",
    "        for key in field_dict:\n",
    "            for ds in field_dict[key]:\n",
    "                if key is not None:\n",
    "                    ds_name=beam+'/land_ice_segments/'+key+'/'+ds\n",
    "                else:\n",
    "                    ds_name=beam+'/land_ice_segments/'+ds\n",
    "                if index is not None:\n",
    "                    D[ds]=np.array(h5f[ds_name][index])\n",
    "                else:\n",
    "                    D[ds]=np.array(h5f[ds_name])\n",
    "                if '_FillValue' in h5f[ds_name].attrs:\n",
    "                    bad_vals=D[ds]==h5f[ds_name].attrs['_FillValue']\n",
    "                    D[ds]=D[ds].astype(float)\n",
    "                    D[ds][bad_vals]=np.NaN\n",
    "        D['data_start_utc'] = h5f['/ancillary_data/data_start_utc'][:]\n",
    "        D['delta_time'] = h5f['/' + beam + '/land_ice_segments/delta_time'][:]\n",
    "        D['segment_id'] = h5f['/' + beam + '/land_ice_segments/segment_id'][:]\n",
    "    if epsg is not None:\n",
    "        xy=np.array(pyproj.proj.Proj(epsg)(D['longitude'], D['latitude']))\n",
    "        D['x']=xy[0,:].reshape(D['latitude'].shape)\n",
    "        D['y']=xy[1,:].reshape(D['latitude'].shape)\n",
    "    temp=file_re.search(filename)\n",
    "    D['rgt']=int(temp['rgt'])\n",
    "    D['cycle']=int(temp['cycle'])\n",
    "    D['beam']=beam\n",
    "    return D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Next we will read in the files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 4 data files of which 0 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "# ATL06_files=glob.glob(os.path.join(datapath, 'PIG_ATL06', '*.h5'))\n",
    "rgt = '0848'\n",
    "ATL06_files=glob.glob(os.path.join(datapath, '*' + rgt + '*.h5'))\n",
    "\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files[:10]:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 9 data files of which 0 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "# ATL06_files=glob.glob(os.path.join(datapath, 'PIG_ATL06', '*.h5'))\n",
    "rgt = '0537'\n",
    "ATL06_files=glob.glob(os.path.join(datapath, '*' + rgt + '*.h5'))\n",
    "\n",
    "#D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files[:10]:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Then, we will plot the ground tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b30845f62f1a4c28aca64a555acb7db4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-887950., -356950.,  183825.,  561825.]), 'origin': 'lower'}\n",
      "{'cmap': 'gray', 'clim': [14000, 17000], 'extent': array([-887950., -356950.,  183825.,  561825.]), 'origin': 'lower'}\n"
     ]
    }
   ],
   "source": [
    "plt.figure(figsize=[8,8])\n",
    "hax0=plt.gcf().add_subplot(211, aspect='equal')\n",
    "MOA.show(ax=hax0, cmap='gray', clim=[14000, 17000]);\n",
    "hax1=plt.gcf().add_subplot(212, aspect='equal', sharex=hax0, sharey=hax0)\n",
    "MOA.show(ax=hax1, cmap='gray', clim=[14000, 17000]);\n",
    "for fname, Di in D_dict.items():\n",
    "    cycle=Di['cycle']\n",
    "    if cycle <= 2:\n",
    "        ax=hax0\n",
    "    else:\n",
    "        ax=hax1\n",
    "    #print(fname)\n",
    "    #print(f'\\t{rgt}, {cycle}, {region}')\n",
    "    ax.plot(Di['x'], Di['y'])\n",
    "    if True:\n",
    "        try:\n",
    "            if cycle  < 3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}\", clip_on=True)\n",
    "            elif cycle==3:\n",
    "                ax.text(Di['x'][0], Di['y'][0], f\"rgt={Di['rgt']}, cyc={cycle}+\", clip_on=True)\n",
    "        except IndexError:\n",
    "            pass\n",
    "hax0.set_title('cycles 1 and 2');\n",
    "hax1.set_title('cycle 3+');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "read 4 data files of which 0 gave errors\n"
     ]
    }
   ],
   "source": [
    "# find all the files in the directory:\n",
    "# ATL06_files=glob.glob(os.path.join(datapath, 'PIG_ATL06', '*.h5'))\n",
    "rgt = '0848'\n",
    "ATL06_files=glob.glob(os.path.join(datapath, '*' + rgt + '*.h5'))\n",
    "\n",
    "D_dict={}\n",
    "error_count=0\n",
    "for file in ATL06_files[:10]:\n",
    "    try:\n",
    "        D_dict[file]=atl06_to_dict(file, '/gt2l', index=slice(0, -1, 25), epsg=3031)\n",
    "    except KeyError as e:\n",
    "        print(f'file {file} encountered error {e}')\n",
    "        error_count += 1\n",
    "print(f\"read {len(D_dict)} data files of which {error_count} gave errors\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Repeat track elevation profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A revised code to plot the elevations of segment midpoints (h_li):\n",
    "def plot_elevation(D6, ind=None, **kwargs):\n",
    "    \"\"\"\n",
    "    Plot midpoint elevation for each ATL06 segment\n",
    "    \"\"\"\n",
    "    if ind is None:\n",
    "        ind=np.ones_like(D6['h_li'], dtype=bool)\n",
    "    # pull out heights of segment midpoints\n",
    "    h_li = D6['h_li'][ind]\n",
    "    # pull out along track x coordinates of segment midpoints\n",
    "    x_atc = D6['x_atc'][ind]\n",
    "\n",
    "    plt.plot(x_atc, h_li, **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Data Visualization**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "filename=/home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190523195046_08480311_003_01.h5, exception=name 'map_ax' is not defined\n",
      "filename=/home/jovyan/shared/surface_velocity/FIS_ATL06/processed_ATL06_20190822153035_08480411_003_01.h5, exception=name 'map_ax' is not defined\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2aa3a9970cda4c53954988a55dc7bf38",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "D_2l={}\n",
    "D_2r={}\n",
    "\n",
    "# specify the rgt here:\n",
    "rgt=\"0027\"\n",
    "rgt=\"0848\" #Ben's suggestion\n",
    "\n",
    "# iterate over the repeat cycles\n",
    "for cycle in ['03','04','05','06','07']:\n",
    "    for filename in glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5')):\n",
    "        try:\n",
    "            # read the left-beam data\n",
    "            D_2l[filename]=atl06_to_dict(filename,'/gt2l', index=None, epsg=3031)\n",
    "            # read the right-beam data\n",
    "            D_2r[filename]=atl06_to_dict(filename,'/gt2r', index=None, epsg=3031)\n",
    "            # plot the locations in the previous plot\n",
    "            map_ax.plot(D_2r[filename]['x'], D_2r[filename]['y'],'k');  \n",
    "            map_ax.plot(D_2l[filename]['x'], D_2l[filename]['y'],'k');\n",
    "        except Exception as e:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "plt.figure();\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylabel('elevation');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Now, we need to pull out a segment and cross correlate:\n",
    "\n",
    "Let's try 2.93e7 through x_atc=2.935e7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9fd4ec2c650f49b8a6e47a9de26d9e81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cycles = [] # names of cycles with data\n",
    "for filename, Di in D_2l.items():\n",
    "    cycles += [str(Di['cycle']).zfill(2)]\n",
    "cycles.sort()\n",
    "    \n",
    "# x1 = 2.93e7\n",
    "# x2 = 2.935e7\n",
    "\n",
    "beams = ['gt1l','gt1r','gt2l','gt2r','gt3l','gt3r']\n",
    "\n",
    "# try and smooth without filling nans\n",
    "smoothing_window_size = int(np.round(60 / dx)) # meters / dx; odd multiples of 20 only! it will break\n",
    "filt = np.ones(smoothing_window_size)\n",
    "smoothed = True\n",
    "\n",
    "### extract and plot data from all available cycles\n",
    "fig, axs = plt.subplots(4,1)\n",
    "x_atc = {}\n",
    "h_li_raw = {}\n",
    "h_li = {}\n",
    "h_li_diff = {}\n",
    "times = {}\n",
    "for cycle in cycles:\n",
    "    # find Di that matches cycle:\n",
    "    Di = {}\n",
    "    x_atc[cycle] = {}\n",
    "    h_li_raw[cycle] = {}\n",
    "    h_li[cycle] = {}\n",
    "    h_li_diff[cycle] = {}\n",
    "    times[cycle] = {}\n",
    "\n",
    "    filenames = glob.glob(os.path.join(datapath, f'*ATL06_*_{rgt}{cycle}*_003*.h5'))\n",
    "    for filename in filenames:\n",
    "        try:\n",
    "            for beam in beams:\n",
    "                Di[filename]=atl06_to_dict(filename,'/'+ beam, index=None, epsg=3031)\n",
    "\n",
    "                times[cycle][beam] = Di[filename]['data_start_utc']\n",
    "                \n",
    "                # extract h_li and x_atc for that section\n",
    "                x_atc_tmp = Di[filename]['x_atc']\n",
    "                h_li_tmp = Di[filename]['h_li']#[ixs]\n",
    "                \n",
    "                # segment ids:\n",
    "                seg_ids = Di[filename]['segment_id']\n",
    "#                 print(len(seg_ids), len(x_atc_tmp))\n",
    "                \n",
    "                # make a monotonically increasing x vector\n",
    "                # assumes dx = 20 exactly, so be carefull referencing back\n",
    "                ind = seg_ids - np.nanmin(seg_ids) # indices starting at zero, using the segment_id field, so any skipped segment will be kept in correct location\n",
    "                x_full = np.arange(np.max(ind)+1) * 20 + x_atc_tmp[0]\n",
    "                h_full = np.zeros(np.max(ind)+1) + np.NaN\n",
    "                h_full[ind] = h_li_tmp\n",
    "                \n",
    "                \n",
    "                x_atc[cycle][beam] = x_full\n",
    "                h_li_raw[cycle][beam] = h_full\n",
    "                              \n",
    "                # running average smoother /filter\n",
    "                if smoothed == True:\n",
    "                    h_smoothed = (1/smoothing_window_size) * np.convolve(filt, h_full, mode=\"same\")\n",
    "                    #h_smoothed = h_smoothed[int(np.floor(smoothing_window_size/2)):int(-np.floor(smoothing_window_size/2))] # cut off ends\n",
    "                    h_li[cycle][beam] = h_smoothed\n",
    "\n",
    "    #                 # differentiate that section of data\n",
    "                    h_diff = (h_smoothed[1:] - h_smoothed[0:-1]) / (x_full[1:] - x_full[0:-1])\n",
    "                else: \n",
    "                    h_li[cycle][beam] = h_full\n",
    "                    h_diff = (h_full[1:] - h_full[0:-1]) / (x_full[1:] - x_full[0:-1])\n",
    "                    \n",
    "                h_li_diff[cycle][beam] = h_diff\n",
    "\n",
    "                # plot\n",
    "                axs[0].plot(x_full, h_full)\n",
    "                axs[1].plot(x_full[1:], h_diff)\n",
    "#                 axs[2].plot(x_atc_tmp[1:] - x_atc_tmp[:-1])\n",
    "                axs[2].plot(np.isnan(h_full))\n",
    "                axs[3].plot(seg_ids[1:]- seg_ids[:-1])\n",
    "\n",
    "\n",
    "\n",
    "        except:\n",
    "            print(f'filename={filename}, exception={e}')\n",
    "\n",
    "\n",
    "            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Speaker: Grace Barcheck**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "64c57d07e4ee4fec8b67ffc2e39c281b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c50f94ddc86b45d6b76e77b17087c62d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1b017af5389c441995b6f11d4418bf68",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "20fdf49a45a74ddca257f8a83a1d3d9d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ebbf00e62d154f19a9dbf2a6453e4bda",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2ce3e0fb6c804f8aac5a4f7de46d21c3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_veloc = len(cycles) - 1\n",
    "\n",
    "segment_length = 3000 # m\n",
    "x1 = 2.935e7# 2.925e7#x_atc[cycles[0]][beams[0]][1000] <-- the very first x value in a file; doesn't work, I think b/c nans # 2.93e7\n",
    "#x1=2.917e7\n",
    "search_width = 800 # m\n",
    "dx = 20 # meters between x_atc points\n",
    "\n",
    "for veloc_number in range(n_veloc):\n",
    "    cycle1 = cycles[veloc_number]\n",
    "    cycle2 = cycles[veloc_number+1]\n",
    "    t1_string = times[cycle1]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t1 = Time(t1_string)\n",
    "                       \n",
    "    t2_string = times[cycle2]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "    t2 = Time(t2_string)\n",
    "    \n",
    "    dt = (t2 - t1).jd # difference in julian days\n",
    "        \n",
    "    velocities = {}     \n",
    "    for beam in beams:\n",
    "        fig1, axs = plt.subplots(4,1)\n",
    "        \n",
    "        # cut out small chunk of data at time t1 (first cycle)\n",
    "        x_full_t1 = x_atc[cycle1][beam]\n",
    "        ix_x1 = np.arange(len(x_full_t1))[x_full_t1 >= x1][0]\n",
    "        ix_x2 = ix_x1 + int(np.round(segment_length/dx))      \n",
    "        x_t1 = x_full_t1[ix_x1:ix_x2]\n",
    "        h_li1 = h_li_diff[cycle1][beam][ix_x1-1:ix_x2-1] # start 1 index earlier because the data are differentiated\n",
    "        \n",
    "        # cut out a wider chunk of data at time t2 (second cycle)\n",
    "        x_full_t2 = x_atc[cycle2][beam]\n",
    "        ix_x3 = ix_x1 - int(np.round(search_width/dx)) # offset on earlier end by # indices in search_width\n",
    "        ix_x4 = ix_x2 + int(np.round(search_width/dx)) # offset on later end by # indices in search_width\n",
    "        x_t2 = x_full_t2[ix_x3:ix_x4]\n",
    "        h_li2 = h_li_diff[cycle2][beam][ix_x3:ix_x4]\n",
    "\n",
    "        # plot data\n",
    "        axs[0].plot(x_t2, h_li2, 'r')\n",
    "        axs[0].plot(x_t1, h_li1, 'k')\n",
    "        axs[0].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        # correlate old with newer data\n",
    "        corr = correlate(h_li1, h_li2, mode = 'valid', method = 'direct') \n",
    "        norm_val = np.sqrt(np.sum(h_li1**2)*np.sum(h_li2**2)) # normalize so values range between 0 and 1\n",
    "        corr = corr / norm_val\n",
    "        \n",
    "        \n",
    "#         lagvec = np.arange( -(len(h_li1) - 1), len(h_li2), 1)# for mode = 'full'\n",
    "#         lagvec = np.arange( -int(search_width/dx) - 1, int(search_width/dx) +1, 1) # for mode = 'valid'\n",
    "        lagvec = np.arange(- int(np.round(search_width/dx)), int(search_width/dx) +1,1)# for mode = 'valid'\n",
    "\n",
    "        shift_vec = lagvec * dx\n",
    "\n",
    "        ix_peak = np.arange(len(corr))[corr == np.nanmax(corr)][0]\n",
    "        best_lag = lagvec[ix_peak]\n",
    "        best_shift = shift_vec[ix_peak]\n",
    "        velocities[beam] = best_shift/(dt/365)\n",
    "\n",
    "        axs[1].plot(lagvec,corr)\n",
    "        axs[1].plot(lagvec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[1].set_xlabel('lag (samples)')\n",
    "\n",
    "        axs[2].plot(shift_vec,corr)\n",
    "        axs[2].plot(shift_vec[ix_peak],corr[ix_peak], 'r*')\n",
    "        axs[2].set_xlabel('shift (m)')\n",
    "\n",
    "        # plot shifted data\n",
    "        axs[3].plot(x_t2, h_li2, 'r')\n",
    "        axs[3].plot(x_t1 - best_shift, h_li1, 'k')\n",
    "        axs[3].set_xlabel('x_atc (m)')\n",
    "        \n",
    "        axs[0].text(x_t2[100], 0.6*np.nanmax(h_li2), beam)\n",
    "        axs[1].text(lagvec[5], 0.6*np.nanmax(corr), 'best lag: ' + str(best_lag) + '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.6*np.nanmax(corr), 'best shift: ' + str(best_shift) + ' m'+ '; corr val: ' + str(np.round(corr[ix_peak],3)))\n",
    "        axs[2].text(shift_vec[5], 0.3*np.nanmax(corr), 'veloc of ' + str(np.round(best_shift/(dt/365),1)) + ' m/yr')\n",
    "\n",
    "        plt.tight_layout()    \n",
    "    fig1.suptitle('black = older cycle data, red = newer cycle data to search across')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_veloc = len(cycles) - 1\n",
    "\n",
    "segment_length = 2000 # m\n",
    "\n",
    "search_width = 800 # m\n",
    "dx = 20 # meters between x_atc points\n",
    "\n",
    "correlation_threshold = 0.65\n",
    "\n",
    "x1 = 2.915e7#x_atc[cycles[0]][beams[0]][1000] <-- the very first x value in a file; doesn't work, I think b/c nans # 2.93e7\n",
    "x1s = x_atc[cycles[veloc_number]][beams[0]][search_width:-segment_length-2*search_width:10]\n",
    "velocities = {}     \n",
    "correlations = {}     \n",
    "\n",
    "for beam in beams:\n",
    "    velocities[beam] = np.empty_like(x1s)\n",
    "    correlations[beam] = np.empty_like(x1s)\n",
    "for xi,x1 in enumerate(x1s):\n",
    "    for veloc_number in range(n_veloc):\n",
    "        cycle1 = cycles[veloc_number]\n",
    "        cycle2 = cycles[veloc_number+1]\n",
    "        t1_string = times[cycle1]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "        t1 = Time(t1_string)\n",
    "\n",
    "        t2_string = times[cycle2]['gt1l'][0].astype(str) #figure out later if just picking hte first one it ok\n",
    "        t2 = Time(t2_string)\n",
    "\n",
    "        dt = (t2 - t1).jd # difference in julian days\n",
    "\n",
    "        for beam in beams:\n",
    "            # cut out small chunk of data at time t1 (first cycle)\n",
    "            x_full_t1 = x_atc[cycle1][beam]\n",
    "            ix_x1 = np.arange(len(x_full_t1))[x_full_t1 >= x1][0]\n",
    "            ix_x2 = ix_x1 + int(np.round(segment_length/dx))      \n",
    "            x_t1 = x_full_t1[ix_x1:ix_x2]\n",
    "            h_li1 = h_li_diff[cycle1][beam][ix_x1-1:ix_x2-1] # start 1 index earlier because the data are differentiated\n",
    "\n",
    "            # cut out a wider chunk of data at time t2 (second cycle)\n",
    "            x_full_t2 = x_atc[cycle2][beam]\n",
    "            ix_x3 = ix_x1 - int(np.round(search_width/dx)) # offset on earlier end by # indices in search_width\n",
    "            ix_x4 = ix_x2 + int(np.round(search_width/dx)) # offset on later end by # indices in search_width\n",
    "            x_t2 = x_full_t2[ix_x3:ix_x4]\n",
    "            h_li2 = h_li_diff[cycle2][beam][ix_x3:ix_x4]\n",
    "\n",
    "            # correlate old with newer data\n",
    "            corr = correlate(h_li1, h_li2, mode = 'valid', method = 'direct') \n",
    "            norm_val = np.sqrt(np.sum(h_li1**2)*np.sum(h_li2**2)) # normalize so values range between 0 and 1\n",
    "            corr = corr / norm_val\n",
    "\n",
    "\n",
    "    #         lagvec = np.arange( -(len(h_li1) - 1), len(h_li2), 1)# for mode = 'full'\n",
    "    #         lagvec = np.arange( -int(search_width/dx) - 1, int(search_width/dx) +1, 1) # for mode = 'valid'\n",
    "            lagvec = np.arange(- int(np.round(search_width/dx)), int(search_width/dx) +1,1)# for mode = 'valid'\n",
    "\n",
    "            shift_vec = lagvec * dx\n",
    "            \n",
    "            if all(np.isnan(corr)):\n",
    "                velocities[beam][xi] = np.nan\n",
    "                correlations[beam][xi] = np.nan\n",
    "\n",
    "            else:\n",
    "                correlation_value = np.nanmax(corr)\n",
    "                if correlation_value >= correlation_threshold:\n",
    "                    ix_peak = np.arange(len(corr))[corr == correlation_value][0]\n",
    "                    best_lag = lagvec[ix_peak]\n",
    "                    best_shift = shift_vec[ix_peak]\n",
    "                    velocities[beam][xi] = best_shift/(dt/365)\n",
    "                    correlations[beam][xi] = correlation_value\n",
    "\n",
    "                else:\n",
    "                    velocities[beam][xi] = np.nan\n",
    "                    correlations[beam][xi] = correlation_value\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4f40c32824e9418095ce272e7d098617",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Along track velocity: all beams')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "ax1 = plt.subplot(211)\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "plt.legend()\n",
    "plt.ylabel('elevation');\n",
    "\n",
    "ax2 = plt.subplot(212,sharex=ax1)\n",
    "for beam in beams:\n",
    "    plt.plot(x1s+dx*(segment_length/2),velocities[beam],'.',alpha=0.2,ms=3,label=beam)\n",
    "plt.ylabel('velocity (m/yr)')\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylim(0,1500)\n",
    "plt.legend()\n",
    "plt.suptitle('Along track velocity: all beams')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### **Median velocity for all 6 beams:**\n",
    "\n",
    "**Above a cross-correlation threshold of 0.65**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a4c3cac0e8d74676955e275bce070aee",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1113: RuntimeWarning: Mean of empty slice\n",
      "  return np.nanmean(a, axis, out=out, keepdims=keepdims)\n",
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1667: RuntimeWarning: Degrees of freedom <= 0 for slice.\n",
      "  keepdims=keepdims)\n",
      "No handles with labels found to put in legend.\n",
      "/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel_launcher.py:43: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1d496ad1d4d84306979c74181753cff9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0.98, 'Correlation values, all beams')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "ax1 = plt.subplot(211)\n",
    "for filename, Di in D_2l.items():\n",
    "    #Plot only points that have ATL06_quality_summary==0 (good points)\n",
    "    hl=plot_elevation(Di, ind=Di['atl06_quality_summary']==0, label=f\"cycle={Di['cycle']}\")\n",
    "    #hl=plt.plot(Di['x_atc'][Di['atl06_quality_summary']==0], Di['h_li'][Di['atl06_quality_summary']==0], '.', label=f\"cycle={Di['cycle']}\")\n",
    "plt.legend()\n",
    "plt.ylabel('elevation');\n",
    "\n",
    "ax2 = plt.subplot(212,sharex=ax1)\n",
    "\n",
    "medians = np.empty(len(x1s))\n",
    "stds = np.empty(len(x1s))\n",
    "\n",
    "for xi, x1 in enumerate(x1s):\n",
    "    corr_vals = []\n",
    "    velocs = []\n",
    "    for beam in beams:\n",
    "        corr_vals += [correlations[beam][xi]]\n",
    "        velocs += [velocities[beam][xi]]\n",
    "    n_obs = len(velocs)\n",
    "    if n_obs >0:\n",
    "        corr_mask = np.array(corr_vals) >= correlation_threshold\n",
    "        veloc_mask = np.abs(np.array(velocs)) < 0.67*segment_length # get rid of segments that are nailed against one edge for some reason\n",
    "        mask = corr_mask * veloc_mask\n",
    "        median_veloc = np.nanmedian(np.array(velocs)[mask])\n",
    "        \n",
    "        std_veloc = np.nanstd(np.array(velocs)[mask])\n",
    "        medians[xi] = median_veloc\n",
    "        stds[xi] = std_veloc\n",
    "        ax2.plot([x1,x1], [median_veloc - std_veloc, median_veloc +std_veloc], '-', color= [0.7, 0.7, 0.7])\n",
    "\n",
    "ax2.plot(x1s, medians, 'k.', markersize=2)\n",
    "\n",
    "# for beam in beams:\n",
    "#     plt.plot(x1s+dx*(segment_length/2),velocities[beam],'.',alpha=0.2,ms=3,label=beam)\n",
    "plt.ylabel('velocity (m/yr)')\n",
    "plt.xlabel('x_atc')\n",
    "plt.ylim(0,1500)\n",
    "plt.legend()\n",
    "plt.suptitle('Median along track velocity')\n",
    "\n",
    "plt.figure()\n",
    "ax1 = plt.subplot(211)\n",
    "for beam in beams:\n",
    "    xvals = x1s+dx*(segment_length/2)\n",
    "    corrs = correlations[beam]\n",
    "    ixs = corrs >= correlation_threshold\n",
    "    ax1.plot(xvals[ixs], corrs[ixs],'.',alpha=0.2,ms=3,label=beam)\n",
    "plt.ylabel('correlation values, 0->1')\n",
    "plt.xlabel('x_atc')\n",
    "plt.legend()\n",
    "plt.suptitle('Correlation values > threshold, all beams')\n",
    "\n",
    "ax1 = plt.subplot(212)\n",
    "for beam in beams:\n",
    "    ax1.plot(x1s+dx*(segment_length/2),correlations[beam],'.',alpha=0.2,ms=3,label=beam)\n",
    "plt.ylabel('correlation values, 0->1')\n",
    "plt.xlabel('x_atc')\n",
    "plt.legend()\n",
    "plt.suptitle('Correlation values, all beams')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison between measures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**848**\n",
    "\n",
    "**Speaker: Lynn Kaluzienski**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![LynnFig](Figure_96.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![comparison1](median_along_track1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**537**\n",
    "\n",
    "**Speaker: Joseph Martin**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![537_velocity](0537_measures.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![comparison](median_along_track2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![537_veloc](537_velocity.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Future Work for the Surface Velocity Team:\n",
    "\n",
    "**Speaker: David Polashenski**\n",
    "\n",
    "- Calculating correlation uncertainty\n",
    "- Considering larger, more complex areas\n",
    "- Pending objectives\n",
    "    - Develop methodology to extract Across Track velocities and test efficacy\n",
    "    - Compare ICESat GLAS methodology (Along Track) to ICESat-2 methodology (Across Track)\n",
    "    - Compare the capabilites of ICESat-2 to extract surface ice velocity from ice shelves and ice streams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:notebook] *",
   "language": "python",
   "name": "conda-env-notebook-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
